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Inexact Shift-and-Invert Arnoldi for Toeplitz Matrix Exponential 


Ting-ting Fen^, Gang WiJ^, Yimin Wej^ 

Abstract 

We revisit the shift-and-invert Arnoldi method proposed in [S. Lee, H. Pang, and H. Sun. Shift-invert 
Amoldi approximation to the Toeplitz matrix exponential, SIAM J. Sci. Comput., 32: 774-792, 2010] for 
numerical approximation to the product of Toeplitz matrix exponential with a vector. In this approach, 
one has to solve two large scale Toeplitz linear systems in advance. However, if the desired accuracy is 
high, the cost will be prohibitive. Therefore, it is interesting to investigate how to solve the Toeplitz 
systems inexactly in this method. The contribution of this paper is in three regards. First, we give a 
new stability analysis on the Gohberg-Semencul formula (GSF) and define the GSF condition number of 
a Toeplitz matrix. It is shown that, when the size of the Toeplitz matrix is large, our result is sharper 
than the one given in [M. Gutknecht and M. Hochbruck. The stability of inversion formulas for Toeplitz 
matrices, Linear Algebra AppL, 223/224: 307-324, 1995]. Second, we establish a relation between the 
error of Toeplitz systems and the residual of Toeplitz matrix exponential. We show that if the GSF 
condition number of the Toeplitz matrix is medium sized, then the Toeplitz systems can be solved in a 
low accuracy. Third, based on this relationship, we present a practical stopping criterion for relaxing 
the accuracy of the Toeplitz systems, and propose an inexact shift-and-invert Arnoldi algorithm for the 
Toeplitz matrix exponential problem. Numerical experiments illustrate the numerical behavior of the 
new algorithm, and show the effectiveness of our theoretical results. 

Keywords: Toeplitz matrix. Matrix exponential, Shift-and-invert Arnoldi, Gohberg-Semencul formula 
(GSF), GSF condition number. 


1 Introduction 

Toeplitz matrices occur in a variety of applications in mathematics and engineering such as complex and 
harmonic analysis, statistics, signal and image processing, information theory, numerical analysis, see [an 
120] and the references therein. In this paper, we are interested in numerical approximation to the product 
of Toeplitz matrix exponential with a vector 


y{t) = exp(-tA)v, (I.l) 

where t is a scalar, v is a given vector, and —tA is a real n x n large Toeplitz matrix whose spectrum is 
located in the left half plane. This problem plays an important role in various application fields such as 
computational finance [13113, numerical solution of Volterra-Wiener-Hopf equations [I], calculating the 
Wiener-Hopf integral equations m, and so on. 
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The Krylov subspace method is an efficient approach to approximate the matrix exponential with a 
vector, especially when the matrix is very large [I1IIH1I23]. Indeed, it is the twentieth dubious way to 
compute the matrix exponential m- In this type of method, the matrix is first projected into a much 
smaller subspace, then the exponential is applied to the projected matrix, and finally the approximation is 
projected back to the original large space [13111123]. This procedure can be achieved by using the Lanczos 
process for symmetric matrices or by the Arnoldi process for non-symmetric matrices, while both procedures 
require only matrix-vector multiplications. 

The shift-and-invert Arnoldi and Lanczos methods were widely investigated to speed up the Arnoldi and 
the Lanczos methods for matrix exponential [Tl[2l. Recently, by making use of the shift-and-invert Arnoldi 
technique, Toeplitz structure and the famous Gohberg-Semencul formula (GSF) d. Pang et al. proposed a 
shift-and-invert Arnoldi method for Toeplitz matrix exponential mui]. An advantage of this approach is 
that it is unnecessary to explicitly form or store the Toeplitz matrix and its inverse, and each matrix-vector 
product can be realized in several Fast Fourier Transformations (FFTs) [H [21]. In the first step of this 
approach, one has to solve two large scale (non-Hermitian) Toeplitz linear systems in a desired accuracy. 
However, if the desired accuracy is very high, the cost for solving the Toeplitz linear systems will be very 
large, especially for some ill-conditioned problems. Thus, it is interesting to investigate how to solve the 
Toeplitz linear systems inexactly in the shift-and-invert Arnoldi method for matrix exponential. 

In this paper, we first give a new stability analysis on the Gohberg-Semencul formula in terms of I-norm 
and 2-norm, and define the “GSF condition number” of a Toeplitz matrix. It is shown that our results are 
sharper than the one given in |12| when the Toeplitz matrix is large. We then establish a relation between 
the error of Toeplitz systems and the residual of Toeplitz matrix exponential. Based on the relationship, we 
present a practical stopping criterion for solving the Toeplitz systems inexactly. 

This paper is organized as follows. In Section 2, we briefly introduce the shift-and-invert Arnoldi method 
for Toeplitz matrix exponential m- In Section 3, we give a stability analysis on the Gohberg-Semencul 
formula and propose an inexact shift-and-invert Arnoldi algorithm. Numerical results given in Section 4 
show the efficiency of our new algorithm and the effectiveness of the theoretical results. 

2 The shift-and-invert Arnoldi method for Toeplitz matrix expo¬ 
nential 

In the shift-and-invert Arnoldi/Lanczos method [3 [13 [13 Hi] |^, the Krylov subspace is constructed by 
using the matrix (/where 7 is a user-prescribed parameter and I is the identity matrix whose order 
is clear from context. Let Vi = v/||v|| 2 , the m-step shift-and-invert Arnoldi process leads to the following 
relation 

(d -f 7 ^) ^rn — -f (^-l) 

where Vm = [vi, V 2 ,..., v^] is an n x m orthonormal matrix, n is the size of the Toeplitz matrix, Hm = 
V^{I + ')A)~^Vm is an m-by-m upper Hessenberg matrix, and is the m-th column of the m-hy-m identity 
matrix. 

Let P = ||v|| 2 , if Hm is invertible, then the shift-and-invert Arnoldi method exploits 
ym(t) = Vm [exp( - (t/y) ■ {Hm - I)) ■ /3ei] = VmVim{t) 
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as an approximation to y{t), where Um{t) = exp( — (t/y) ■ {H^ — /)) • /3ei. The residual is [5] 


rm (0 


- y'^(t) = -AVmUmit) - VmU'^it) 

hm+l,m T TT—1 tj.\ It I \ 

- ^ruHru Um(i) ’ (/ + 7^)Vm+l, 


and 


km(i )||2 = 


^m+l,m T rr—1 

-^I 


At) 


11(7 + 7A)vm+i||2, 


( 2 . 2 ) 


which can be used as a cheap stopping criterion in practice. 

In the TO-step shift-and-invert Arnoldi method, we have to compute m Toeplitz matrix-vector products 
(/ -I- yA)~^Wi, i = 1, 2,..., m. Since 7 is a given shift, we are interested in computing (J -|- 7 ^.)“^ once 
for all. One option is to compute the inverse by some direct methods such as the LU decomposition m- 
However, Toeplitz matrix is often dense, and the computation of the inverse of a large dense matrix is 
prohibitive, especially when the matrix is large. Fortunately, as / -I- yA is also a Toeplitz matrix, we have 
the Gohberg-Semencul formula (GSF) for its inverse. Indeed, the inverse of a Toeplitz matrix T can be 
reconstructed from its first and last columns. More precisely, denote by ei,e„ the first and the last column 
of the n-hy-n identity matrix, and let x = [^O; Ci; • ■ •; Cn-i]"*" and y = [ 770 , 771 ,..., 77„_i]^ be the solutions of 
the following two Toeplitz systems 

Tx = ei and Ty = e„. (2.3) 

If ^0 A Oj then the Gohberg-Semencul formula can be expressed as 
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(2.4) 


where L^, are lower Toeplitz matrices, and Ry, are upper Toeplitz matrices. Consequently, the Toeplitz 
matrix-vector product {I + yA)~^Vi can be realized in several FFTs of length n [THIIII]. We are in a position 
to present the following algorithm for the Toeplitz matrix exponential; for more details, refer to [El- 


Algorithm 1. An shift-and-invert Arnoldi algorithm for product of Toeplitz matrix exponential 
with a vector 

Step 1. Solve the Toeplitz systems (J -I- yA)yi. = ei and (/ -|- 7 A)y = e„; 

Step 2. Choose a convergence tolerance tolexp o,nd the starting vector Vi = v/||v|| 2 ; 

for 7=1,2,... do 

Perform the shift-and-invert Arnoldi process in which the Toeplitz matrix vector products 
(/ -|- yA)~^'Vi are realized through FFTs. 

//||ri(t )||2 < tolexp, then form the approximation yi{t) = ViUi{t) and Stop, else Continue; 

end for 
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In Step 1 of this algorithm, we have to solve two large scale non-Hermitian Toeplitz linear systems (1131). 
If the desired accuracy is too high, then we have to pay a large amount of computational cost for solving 
the Toeplitz linear systems, especially for some ill-conditioned problems. It is interesting to investigate how 
to solve the Toeplitz systems inexactly nani]. 


3 An inexact shift-and-invert Arnold! algorithm for Toeplitz ma¬ 
trix exponential 

In this section, we consider how to solve the Toeplitz systems inexactly in the shift-and-invert Arnoldi 
method. As we solve the Toeplitz linear systems once for all, it can be understood as an “inexact” inverse 
technology. We first give a new stability analysis on the Gohberg-Semencul formula with respect to 1- 
norm and 2-norm, and then establish a relation between the error of Toeplitz systems and the residual 
of Toeplitz matrix exponential. Based on these theoretical results, we propose an inexact shift-and-invert 
Arnoldi algorithm for Toeplitz matrix exponential. 


3.1 A new stability analysis on the Gohberg-Semencul formula and the GSF 
condition number 

In this subsection, we give a stability analysis on the Gohberg-Semencul formula and define the “GSF 
condition number” of a Toeplitz matrix. Let x = [^O)Cij • ■ • >G-i]"*" y = [7jo,rji,... be the 

numerical solutions of Tx = ei and Ty = e„, respectively. If ^ 0, we denote 
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(3.1) 


which is a perturbation to the Toeplitz inverse T The following theorem gives an error analysis on the 
Gohberg-Semencul formula in terms of 1-norm. 

Theorem 3.1. Given e > 0, if ^ 0, 0, let e = ^6 ^^6 relative error of l/^o with respect 

to 1 /^ 0 ) o,nd 

(3.2) 
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[£-b(e-b(l-b£)e)(l-fe)] • min{||x||i, ||y||i}. 
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(3.4) 






























Proof. It follows from (12.41) and (13.11) that 
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Moreover, we have that 
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On the one hand, we obtain from (13.21) that 
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where e = is the relative error of 1/Co- On the other hand, we note from (j2.4p and p.ip that 


and 


From dSU-dSH), we obtain 


||L4li = l|x||i, P,||i = ||y||i<(i + e)lly||i, 
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Similarly, for the second part of (13.51) . we can prove that 

[e + (e + (1 + £)?)(1 + £)] • ||x||i 



-Ll(—Rl) 

< 

1 

Ho ' 

^VCo / 

1 

Co 


(3.10) 


(3.11) 


Combining (13.51) . (I3.10|) and (I3.11L we arrive at 
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By (I2.3L we have that ||x||i < ||r ^||i and ||y||i < ||r ^||i. Thus, 
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a combination of which yields (El)- 

Furthermore, we have the following corollary on the relative error of Toeplitz inverse. 
Corollary 3.1. Under the above notations, there holds 
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□ 

By (13.121) . (|lT||i||y||i)/(|^o|/||x||i) is an enlarge factor of the solution of T~^ over the vector error £. So 
we can give the following definition on the condition number of a Toeplitz matrix. 


Definition 3.1. We define 
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as the 1-norm “GSF condition number” of a Toeplitz matrix. 

Note that |Co|/||x||i is the “proportion” of |^o| with respect to ||x||i, and 
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the “effective” 1-norm condition number of Ty = e„ defined in [Mila, and we have that 

K,{T)=KUT)-K^f,{T), 


(3.13) 


(3.14) 


where ki(T) = ||T||i ||T is the “classical” 1-norm condition number [TT] of the matrix T. 

Remark 3.1. In terms of Corollary \ 3. 11 is an estimation to ki{T). By \3.14\ , ||x||i/|^o| can be 

used as an approximation to Kgg(T). Thus, an advantage of L‘1.13[\ is that one can evaluate the “classical” 
condition number ki{T) of a Toeplitz matrix, and the “effective” condition numbers /vgfj(T), Kgg(T) via 
solving Toeplitz systems, with no need to form the Toeplitz inverse explicitly. 
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Since Vm is orthonormal, it is desirable to investigate the absolute and relative errors of T ^ with respect 
to T~^ according to 2-norm. We have the following result. 


Theorem 3.2. Under the assumptions and notations of Theorem \S. R we have that 
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On the other hand, we can give an upper bound on T~^ — T~^ 

[& -I- (e -I- (1 -I- e)£)(l + £)] • ||x||i||y||i, 
whose proof is similar to that of Theorem 3.1, see [7]. So we have from (I3.17L (13.31) and (I3.18P that 
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Combining (13.151) and (I3.19L we drive 
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a combination of the above two inequalities yields (13.161) . 
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Remark 3.2. In ( \1‘A p.321f), Gutknecht and Hochbruck analyzed the stability of the Gohberg-Semencul 
formula, and gave the following two upper bounds for the absolute and relative errors with respect to T~^: 
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where n is the size of the Toeplitz matrix, e is the machine precision and e is a normwise relative error bound 
that satisfies 

Plla < ||x||2 • (1 + e), ||y||2 < l|y||2 • (1 +£)• 


By setting e = 0, and reduce to 
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respectively. 

Compared with l\3.22^ and the new bounds given in and <\3.15\\ are independent of n, while 

(I.V.iffI) relies on instead of n. Thus, our new upper bounds can be sharper than II .V . and (l,V. 28^ when 
n is large and e = 0{s). See Example 8 of Section 4 for a comparison of these upper bounds. 


3.2 Relationship between the error of Toeplitz systems and the residual of 
Toeplitz matrix exponential 

In this subsection, we establish a relationship between the error of Toeplitz systems and the residual of 
Toeplitz matrix exponential, and propose an inexact shift-and-invert Arnoldi method for product of a Toeplitz 
matrix exponential with a vector. It is shown that if the GSF condition number of the Toeplitz matrix is 
medium sized, we can solve the Toeplitz systems in a relatively low accuracy. 

For simplicity, in the following we denote T = I + yd whenever necessary. Indeed, if the Toeplitz 
systems (ESI) are solved inexactly, the errors of the matrix-vector products can be expressed as = T — 
T“^Vi, i = I, 2,..., m. Let Fm = [fi, f 2 ,..., fm], we get the following relation for the m-step “inexact” 
shift-and-invert Arnoldi procedure 

(/+ yd) + Em = VmHm + hjn+l,m'Vm+lG^, (3.24) 

where Kn = [vi, V 2 ,..., v^] is an n x m orthonormal matrix, and Hm = [(/ -|- yd)“^ + FmV.^] Vm is an 

m X m upper Hessenberg matrix. Note that Vm is different from the one given in (EH), and the subspace 
spanned by Vm is not a Krylov subspace any more. 

Lemma 3.1. If Hm is invertible, denote G = —^{I + l2^)FmHfifVm, then the “inexact” shift-and-invert 
Arnoldi relation \3.24\ can be rewritten as 

(d + G)Vm = Vm (h-^ - ^)) - + yd)v„+ieTif-i. (3.25) 

Proof. Multiplying I -\- yd on both sides of (I3.24|) yields 

Vm + (4 + ^A)Fm = {I + lA)VmHm + hm+l,m{I + 

that is, 

AVm - -(/ + lA)FmH-^ = -Vm (l - Hm) H'^ - + yd)V„+i6^ 

'y 'j \ / ^ 

= -Vm (h-^ - /) - 1^22+hlE^i + yd)v™+ieTiJ-i 
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The above equation can be rewritten as 


+ lA)FmH-^V^^ Vm = {A + G)Kn = Kn (i?”' “ + 7 A)v„+ieTif’i. 

□ 

Let Um{t) = exp( — (t/q) • — /)) • /3ei, then we can use ym{t) = VmUm{t) as an approximation to 

y{t). The “real” residual is defined as [5] 

^real ^ _AV^U^(t) - V^U^t) ■ (3.26) 

However, it is not computable since AVm is unavailable in practice. Thus, we define 
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as the “computed” residual. Moreover, 
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(3.28) 


which can be used as a cheap stopping criterion in the “inexact” shift-and-invert Arnoldi method for Toeplitz 
matrix exponential. 

We are ready to provide a practical stopping criterion for solving the Toeplitz systems inexactly. The 
key is how to investigate the distance between r”®®* and It is seen that 

nreal_^compU = || GKnU„ (t) || 2 
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[e+ (e + (1 + e)e)(l +e)] • ||x||i||y||i, i = 1,2,..., m. 
If £ < e <C 1, then £ + (s + (1 + £)£)(! + £) < 3£ + G(£^), and 
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(3.30) 


where we omit the high order term G(£^). 
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Denote H^\im{t) = [ai,a 2 , ■ • ■ ,ctm]'^, then 




m 




2 


< max ||f ,||2 • 

l<2<m 

2 

Hm^Mm(t) 


(3.31) 


< max \\ii\\ 2 -y/rn 

l<i<m 


Combining (I3.29I1 . (13.301) and (13.3111 . we have that 
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Let to/exp be the convergence threshold for the shift-and-invert Arnoldi method for solving (jl.lll . If 
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In conclusion, we have the following theorem. 

Theorem 3.3. Under the above notations and assumptions, if 

ICo|/||x||i 


e < 



'TI * tolejQp 

6^/m ■ 


2 l|Um(t)||2 


I7I • II/ + 7AII1 • tol 


then 


ICo|/l|x||i _ 

ll^ + 7^lli||y||i 6V^||/ + 7 A|| 2 - 

|, reaJ_j.comp|, , 


exp 


h;, 


|Um(t)||5 


(3.32) 


Remark 3.3. We notice that 


(7 + 771) = 


ll^ + 7^lli||y|| 


ICo|/||x||i 

is just the 1-norm “GSF condition number” defined in Definition \3.1[ which can be utilized as an estimation 
to the 1-norm condition number of I + 7 A. Furthermore, ()7. 3 ‘3 can be reformulated as 


e < 


I 7 I • ll^ + 7^lli -to/e 


kP^(/ + 7^) 6V^||/ + 7A||5 


Hj, 


|Um(t)|h 


(3.33) 


This implies that if the GSF condition number of the Toeplitz matrix is medium sized, we can solve the Toeplitz 
systems in a (relatively) low accuracy. Otherwise, we have to solve the Toeplitz systems in a (relatively) high 
accuracy. 


Remark 3.4. Unfortunately, the parameters 


H. 


-1 


, ||um(t)|| 2 , and (I + ^A) are unavailable a prior. 


We notice that 


H- 


is uniformly bounded and ||um(t )||2 = C)(||y(t)|| 2 ) as the shift-and-invert Arnoldi 
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method converges. Therefore, if + ^A) is medium sized and ||/ + 7 AII 1 = O ^ 


H- 


we 


suggest using 


|r.||2,||r,||2}< 


I 7 I 


• toL 


(3.34) 


6y/m- max{||fcol|| 2 , ||frow|| 2 } 
as the stopping criterion for solving the Toeplitz systems, where r^; = ei — Tx and Vy = Bn — Ty are the 
residuals of the Toeplitz systems, and fcol, frow are the first column and first row of I + yA, respectively. 

In summary, we propose the following “inexact” shift-and-invert Arnold! algorithm for solving the Toeplitz 
matrix exponential problem (HU. 


Algorithm 2. An inexact shift-and-invert Arnoldi algorithm for product of Toeplitz matrix 
exponential with a vector 

This algorithm is similar to Algorithm 1, except for the Toeplitz linear systems are solved inexactly in 
Step 1, with the stopping criterion described in l\3.S4l . 

We point out that several results on computing matrix functions using inexact Krlov methods have been 
developed in other contexts, where 2-norm estimates are given. In [5], Frommer et al. considered how to 
cheaply recover a secondary Lanczos process starting at an arbitrary Lanczos vector. This secondary process 
is then used to efficiently obtain computable error estimates and error bounds for the Lanczos approximations 
to the action of a rational matrix function on a vector, e.g., the matrix sign function. 


4 Numerical experiments 


In this section, we perform some numerical examples to show the efficiency of Algorithm 2 and the effective¬ 
ness of our theoretical results. All the numerical experiments were run on two core Intel(R) Core(TM)2 E7400 
processor with CPU 2.8 GHz and RAM 1.99 GB, under the Windows 7 operating system. The experimental 
results were obtained by using a MATLAB 7.7 implementation with machine precision e ~ 2.22 x 10“^®. 
two core Intel(R) Gore(TM)2 E7400 processor with GPU 2.8 GHz and RAM 2 GB 
As was done in [ID, we use the (unrestarted) GMRES algorithm [21] with T. Chan’s optimal (circulant) 
preconditioner EiiBiin] for solving the Toeplitz systems in Algorithm 1 and Algorithm 2. Let tolexp, tolgys 
be the tolerance for computing the Toeplitz matrix exponential-vector product, and that for solving the 
Toeplitz systems, respectively. Denote by A4 the optimal preconditioner due to T. Chan, by q = x (or y) 
the approximate solution of the Toeplitz system, and by b = ei (or e„) the right-hand side. In Algorithm 2 
we use 


\M-^h - M~^{I + yA)^ < 


ItI 

-;-- • tolexxD = tolsys (4.1) 

evlOO • max{||fcol||2, ||frow||2} 


as the stopping criterion for the Toeplitz systems, where fcol and frow denote the first column and the first 
row of J -|- 7 A, respectively. This algorithm mimics solving the two Toeplitz systems “inexactly” with an 
iterative solver. 

In Algorithm 1, we use 


|7W”ib - M-\l + 7A)q||2 < 10“^^ = tols 


(4.2) 


as the stopping criterion for the Toeplitz systems. This algorithm mimics solving the Toeplitz systems 
“exactly” via an iterative solver. Let y(t) be the “exact” solution and ym{t) be the approximate solutions 
obtained from Algorithm 1 or Algorithm 2, then we define 

\\y{t) - ym{t )\\2 


Error = 


\\yit)h 


(4.3) 
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as the relative error of the approximation ym{t)- Except for Example 1, the “exact” solution y(t) is calculated 
by using the MATLAB bulit-in function expm.m. In the tables below, we denote by CPU the CPU time in 
seconds. We choose the vector v = [1,1,...,!]"'" for all the numerical experiments in this section. 

Example 1. In this example, we aim to show the effectiveness of our inexact strategy (|3.34ll . as well 
as the superiority of Algorithm 2 over Algorithm 1. The Toeplitz matrix A is generated by the even 
function 0^ defined on [— tt, tt]. We want to compute y(t) = exp(—tA)v with t = l,tolexp = 10“® and 
n = 1 X 10^,2 X 10®,... ,5 X 10®, respectively. Since the size n of the Toeplitz matrix is very large, the 
MATLAB build-in function expm.m is infeasible for this problem. As a compromise, we run Algorithm 1 
with the convergence tolerance to^exp = 10“^"'^ for the “exact” solution y(t). Table 1 lists the numerical 
results. 

We see from Table 1 that Algorithm 2 converges much faster than Algorithm 1 in practical calculations, 
and the inexact strategy is both efficient and reliable. Thanks to (13.3411 . it is only necessary to solve the 
Toeplitz system in the accuracy of 0(10“®) instead of 10“^^. Furthermore, the approximate solutions 
computed from the two methods have the same accuracy in terms of Error. 


n 

Algorithm 

tolgys 

Error 

CPU 

1 X 10® 

Algorithm 1 

1.000 X 10“*'* 

4.615 X 10“*' 

2.531 


Algorithm 2 

1.239 X 10“® 

4.615 X 10“*' 

1.297 

2 X 10® 

Algorithm 1 

1.000 X 10“*4 

3.263 X 10“^ 

7.360 


Algorithm 2 

1.239 X 10“® 

3.263 X 10“^ 

3.422 

3 X 10® 

Algorithm 1 

1.000 X 10“*'* 

2.664 X 10“*' 

13.375 


Algorithm 2 

1.239 X 10“® 

2.664 X 10“*' 

6.031 

4 X 10® 

Algorithm 1 

1.000 X 10“*4 

2.307 X 10“*' 

20.313 


Algorithm 2 

1.239 X 10“® 

2.307 X 10“*' 

9.125 

5 X 10® 

Algorithm 1 

1.000 X 10“*4 

2.064 X 10“^ 

27.719 


Algorithm 2 

1.239 X 10“® 

2.064 X 10“*' 

10.187 


Table 1, Example 1: A comparison of Algorithm 1 and Algorithm 2, t = 1, 7 = 1/10 and tolexp ~ 10 ®. 

Example 2. The aim of this example is two-fold. First, we show the effectiveness of Theorem 13.31 
Second, we illustrate that our proposed 1-norm “GSF condition number” ()3.13ll is a good estimation to the 
1-norm “classical condition number” of a Toeplitz matrix. For the first aim, we run Algorithm 2 with the 
stopping criterion tolexp chosen as 10 “®, 10 “'*,..., 10 “*®, and try to show that 

||r”“'-r™™P||2=(!l(toUp). 

In order to compute the “real” residual, we first form the approximation ym(t) = UmUm(t) explicitly, and 
then compute by (13.261) . The convergence threshold tolsys for the Toeplitz systems is determined by 
using (13.341) . 

There are two test problems in this example, both of which are from |15j . The first test matrix is the 
non-Hermitian Toeplitz matrix generated by the function f{9) = 9^ + i- 9^,i = y/-A, 9 G [—tt,tt]. Notice 
that Re(f) = d® > 0 is an even function, and Im(f) = 0® is an odd function. Table 2 lists the numerical 
results of Algorithm 2 for exp(— tA)v with t = 1, 7 = I/IO, and n = 3000. For the first test problem, 
we have + 7 A) Ri 1.275 x 10®, which is of medium sized. In the second test problem, we consider 

pricing options for a single underlying asset in Merton’s jump-diffusion model [lain]. As the real part 
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of the eigenvalues of the Toeplitz matrix are less equal to zero, we are interested in computing exp(fA)v 
with t > 0, for more details, see Example 3 of [15]. Table 3 gives the numerical results of Algorithm 2 with 
t = 1, 7 = 1, and n = 3000. For this test problem, we have + 7 A) Ri 6.296 x lO"^, which is relatively 

large. 

Two remarks are in order. First, we see that — r ™™^||2 and tole^p are about in the same order in 

all the cases. This illustrates the effectiveness of Theorem 13.31 as well as the efficiency of the inexact strategy 
(I3.34p . Second, we observe from Table 2 and Table 3 that, if the GSF condition number {I + 7 A) is 
medium sized, one can solve the Toeplitz systems ( 1 ^ with a relatively low accuracy. Otherwise, we have 
to solve them with a relatively high accuracy. For instance, if we choose toZexp = 10“®, one has to solve 
the Toeplitz systems in an accuracy of 0(10“®) for the first test problem, while an accuracy of 0(10“^^) is 
required for the second test problem. 


tolgxp 

tolgys 


Error 

10“® 

1.010 X 10-6 

1.519 X 10-3 

2.679 X 10-4 

10-4 

1.010 X 10“^ 

5.352 X 10-6 

6.480 X 10-6 

10-6 

1.010 X 10“® 

4.839 X 10-6 

3.985 X 10-6 

10“® 

1.010 X 10-44 

1.679 X 10-3 

1.701 X 10“® 

10-10 

1.010 X 10-43 

2.667 X 10“® 

2.607 X 10-4® 


Table 2, the 1st test problem of Example 2: Numerical results of Algorithm 2 with different tolexp for 
computing e“*^v, t = 1, 7 = 1/10, n = 3000; ftf+ 7 A) « 1.275 x 10^. 


t'Olexp 

tolgys 


Error 

10“® 

4.236 X 10“® 

6.958 X 10-® 

9.375 X 10-6 

10-4 

4.236 X 10-44 

7.347 X 10-4 

6.817 X 10-'4 

10-6 

4.236 X 10-43 

3.320 X 10-6 

3.056 X 10"® 

10-3 

4.236 X 10-46 

1.417 X lO-^" 

2.364 X 10-44 

0 

7 

0 

1—1 

4.236 X 10-44’ 

5.105 X 10-44 

2.021 X 10-44 


Table 3, the 2nd test problem of Example 2: Numerical results of Algorithm 2 with different tolexp for 
computing e*^v, t = 1, 7 = 1, n = 3000; + 7 A) « 6.296 x lO’^. 


When n = 1000, 2000,3000 and 4000, we list in Table 4 the 1-norm GSF condition number -|- 7 A)( 

where we use max ||fcol|| 1, ||frow||i instead of ||/ -I-7AII 1 ), the 1-norm classical condition number ki{I -I-7A) 
(evaluated by using the MATLAB command cond{I -h 7A, 1)) and its estimation -|- 7A) (evaluated 

by using the MATLAB command condest{I + 7A)); as well as the GPU time in seconds for solving them 
(in brackets). It is seen that the GSF condition number is about one to two times larger than the classical 
condition number, and the former is a good estimation to the latter. Furthermore, the CPU time for 
^GSF(j is much less than that for ki(/ -I-7A) and -I-7A), especially when n is large. Thus, the 

1-norm GSF condition number is a competitive alternative to the classical condition number for Toeplitz 
matrices. 

Example 3 . In this example, we try to show that our new bounds ()3.I5|) and ()3.I6I1 are sharper than 
(13.221) and (13.231) . The test matrix is the “gallery” matrix generated by the MATLAB command A = 


13 


















Test problem 

n 

Ki{I -I-7A) 

Kf^I + jA) 

+ -fA) 


1000 

65.284(0.219) 

65.284(0.109) 

79.037(0.031) 

1st test problem 

2000 

89.546(1.469) 

89.546(0.562) 

1.071 X 102(0.141) 


3000 

1.073 X 102(4.641) 

1.073 X 102(1.547) 

1.275 X 102(0.265) 


4000 

1.218 X 102(10.594) 

1.218 X 102(3.328) 

1.442 X 102(0.500) 


1000 

2.436 X 10^(0.219) 

2.436 X 10®(0.094) 

6.989 X 10®(0.0620) 

2nd test problem 

2000 

9.736 X 10^(1.469) 

9.736 X 10®(0.531) 

2.797 X 10^(0.218) 


3000 

2.190 X 10^(4.672) 

2.190 X 10^(1.531) 

6.296 X 10^(0.453) 


4000 

3.893 X 10^(10.656) 

3.893 X 10^(3.344) 

1.119 X 10®(0.828) 


Table 4, Example 2 : The values of + 'fA), + 7^4 ), ki(7 + 77 I) and the CPU time in seconds for 

computing them (in brackets), n = 1000, 2000,3000,4000. 


galleryCparter',n) |14] . It is a Toeplitz matrix whose singular values are close to tt. Let x and y be the 
“exact” solutions of the systems (7 + 77l)x = ei and (7 + 7A)y = e„, respectively, which are computed from 
running the preconditioned (unrestarted) GMRES algorithm with tolgys = Then we form x in the 

following way: 

f = randn(n, 1); f = f/||f|l; x = x + £:||x|| • f; 

where randn{n, 1) is a vector of length n with normally distributed random entries, and || • || is 1-norm 
for (13.151) and (13.161) . and 2-norm for (13.221) and (I3.23p . The vector y is formed in a similar way. In this 
example, we choose e = 10“®, 10“®, 10“^^ and n = 1000,2000,3000,4000, respectively. In order to show 
the sharpness of our results, we also present the “exact” absolute and relative errors ||T“^ — and 

—prrT|i ——■ Tables 4 and 5 report the numerical results. It is seen that our upper bounds are sharper than 
those due to Gutknecht and Hochbruck, especially when n is large. 


e 

n 




rp — 1 _ ^—1 


(13.151) 

(13.221) 


1000 

1.114 X 10-*i 

3.358 X 10-5 


3.238 X 10- 

6 

10"® 

2000 

1.217 X 10-*i 

6.717 X 10-5 


3.672 X 10- 

6 


3000 

1.281 X 10-5 

1.007 X 10-2 


3.948 X 10- 

6 


4000 

1.329 X 10-5 

1.343 X 10-2 


3.767 X 10- 

6 


1000 

1.113 X 10-® 

3.358 X 10-6 


3.182 X 10- 

9 

10-9 

2000 

1.218 X 10-® 

6.717 X 10-6 


3.414 X 10- 

9 


3000 

1.282 X 10-® 

1.007 X 10-5 


3.947 X 10- 

9 


4000 

1.328 X 10-® 

1.343 X 10-5 


3.961 X 10- 

9 


1000 

1.113 X 10-11 

3.358 X 10-9 

3.385 X 10-12 

10-12 

2000 

1.217 X 10-11 

6.717 X 10-9 

4.102 X 10-12 


3000 

1.282 X 10-11 

1.007 X 10-5 

3.930 X 10-12 


4000 

1.328 X 10-11 

1.343 X 10-5 

4.071 X 10-12 


Table 5, Example 3: A comparison of the absolute error bounds (13.151) and (13.221) . t = 1,^ = 1/10, 
e = 10“®, 10"^, 10"^^ and n = 1000, 2000, 3000, 4000. 
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£ 

n 

(13.161) 

(13.231) 

l|T-i|U 


1000 

3.524 X 10-4 

4.813 X 10-° 

2.932 X 10-° 

10-® 

2000 

5.447 X 10-4 

9.642 X 10-° 

3.739 X 10-° 


3000 

7.018 X 10-4 

1.447 X 10-2 

3.829 X 10-° 


4000 

8.403 X 10-4 

1.931 X 10-2 

3.652 X 10-° 


1000 

3.519 X 10-^ 

4.813 X 10-° 

3.626 X 10-° 

10-® 

2000 

5.443 X 10-^ 

9.642 X 10-° 

3.328 X 10-° 


3000 

7.018 X 10-^ 

1.447 X 10-° 

3.618 X 10-° 


4000 

8.403 X 10-^ 

1.931 X 10-° 

3.718 X 10-° 


1000 

3.521 X 10-4° 

4.813 X 10-° 

3.427 X 10-42 

10-12 

2000 

5.443 X 10-4° 

9.642 X 10-° 

3.964 X 10-42 


3000 

7.017 X 10-4° 

1.447 X 10-8 

4.893 X 10-42 


4000 

8.402 X 10-4° 

1.931 X 10-8 

3.991 X 10-42 


Table 6, Example 3: A comparison of the relative error bounds (13.1611 and (13.231) . t = 1,7 = 1/10, 
e = 10“®, 10"'^, 10"^^ and n = 1000, 2000, 3000, 4000. 


5 Conclusion 

In this paper, we analyze and further develop an inexact shift-and-invert Arnoldi method for the problem 
of numerical approximation to the product of Toeplitz matrix exponential with a vector. First, we give an 
improved stability analysis on the Gohberg-Semencul formula (GSF) for the inverse of a Toeplitz matrix, 
and our result is independent of the size of the matrix in question. Moreover, we define the “GSF condition 
number” of a Toeplitz matrix. An advantage is that we can evaluate the “classical” condition number and 
the effective condition numbers of a Toeplitz matrix via solving Toeplitz systems, with no need to form 
the Toeplitz inverse explicitly. Second, we establish a relation between the error in approximating Toeplitz 
systems and the residual of its matrix exponential. Third, we provide a practical stopping criterion for the 
accuracy in approximating the Toeplitz systems in the inexact shift-and-invert Arnoldi algorithm for Toeplitz 
matrix exponential. It is shown that if the I-norm “GSF condition number” -I-7A) is medium sized, 

then the Toeplitz systems can be solved in a relatively low accuracy. 
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